arXiv:1504.07908vl [cs.PF] 29 Apr 2015 


Inhomogeneous CTMC Model of a Call Center 
with Balking and Abandonment 

Maciej Rafal Burak 
Applied Informatics, 

West Pomeranian University of Technology 
ul. Sikorskiego 37, Szczecin, Poland 

http: //www. we. zut. edu. pi 

April 30, 2015 


Abstract 

This paper considers a nonstationary multiserver queuing model with 
abandonment and balking for inbound call centers. We present a contin¬ 
uous time Markov chain (CTMC) model which captures the important 
characteristics of an inbound call center and obtain a numerical solu¬ 
tion for its transient state probabilities using uniformization method with 
steady-state detection. 

Keywords: call center, transient, Markov processes, numerical meth¬ 
ods, uniformization, abandonment, balking 


1 Introduction 


The problem of managing operations of a telephone call center in an efficient 
way has a long history in the area of operational research and is a topic of 
current research in various disciplines (see e.g. Aksin et al. | |200^ 


or 


Cans et al. 


2003 for extensive overviews). From the modeling point of view they can be 


viewed as queuing systems. 

Such a queuing model can be described by a corresponding continuous time 
Markov chain (CTMC) whose steady-state distribution can be easily deter¬ 
mined, either analytically - with the Erlang-C formula for the simplest M/M/n 
model or with the Erlang-A formula for its version augmented with exponen¬ 
tial patience time as proposed in Brown et al. [2005 ; or numerically for more 


complicated models (as in Deslauriers et al. 2007 or recently Phung-Duc and 
Kawanishi |2014| ). However, as real call centers are time inhomogenous, with 


varying arrival rates and changing number of servers - scheduled to meet the 
forecasted demand and in order to provide break time, stationary models cannot 
be applied directly. It is, therefore, common to use approximations, assuming 
the system being pointwise stationary. Examples of such well established meth¬ 
ods can be found e.g. 


Green et al. 

|2007|, 

Aksin et al. |2007 or in 

Brown 


et al. [2005 . Unfortunately, stationary approximations are in many cases not 


adequate. For example, Deslauriers et al. [2007 compared them with simula¬ 


tions based on real inbound call center data, with the conclusion that due to 
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the nonstationarity only some of the performance measures can be estimated 
with satisfactory accuracy. Ingolfsson in |Ingolfsson et ah |2010 compared them 
with an inherently transient model and found their results significantly inaccu¬ 
rate or even entirely unreliable. Despite this, their widespread use is commonly 
justified by simple implementation and low computational costs. 

Many authors proposed to use simulation, which can achieve any desired 
accuracy. However, in order to achieve acceptable precision, very long com¬ 
putational times are needed, which makes it often impracticable for common 
applications like schedule planning. 

An alternative approach, which is very effective in terms of the accuracy 
of the model, is to analyze transient CTMC using numerical methods, solving 
effectively their corresponding system of ordinary differential equations (ODE’s) 
as proposed in Ingolfsson et al. (2007 , Bylina et al. (2009 or by the author in 
Burak] (MI] 


Other, less computationally intensive, analytical methods that can approx¬ 
imate such nonstationary systems more accurately than stationary models are 
closure approximations and fluid and diffusion approximations, discussed e.g. in 


Brown et al. 

[2005 , 

Green et al. 

[2(1(17 . 

Czachorski et al. 

[2009 

et al. [2014 

or, for the direct comparison of some examples c 


with the numerical methods and stationary approximations, in [Ingolfsson et ah] 


2007 


Although there is a number of papers dealing with the phenomena of cus¬ 
tomer balking and abandonment in multiserver queues (e.g. Brown et al. [2005| |^ 
delbaum and Zeltyn [2009 Whitt [2006 Artalejo and Pla [2009 or recently 
Phung-Duc and Kawanishi [2014[), they concentrate on stationary models or 


approximations. To the best of our knowledge, an inherently transient CTMC 
model dealing with both balking and abandonment of a call center, has never 
been investigated. 

The main objective of this work is to model such non-stationary systems, 
using transient analysis of corresponding CTMC, in a reliable and precise way, 
with computational efficiency enabling its use for practical applications - in 
particular, as a much more accurate replacement to the Erlang-C and Erlang-A 
formulas, used by practitioners for quantitative call center management. 

In this paper we model an inbound telephone call center with balking and 
abandonment, i.e. the customer may not stay in the queue once realizing he is 
put on hold, or abandon the queue if the waiting time is too long, extending the 


nonstationary M/M/n queuing model analyzed by the author in Burak 2014 . 

The paper is structured as follows. In the next section the model and the 
basic notation are introduced. Section 3 reviews the proposed multi-step uni- 
formization algorithm with steady-state detection and section 4 presents the 
results of numerical experiments. The paper ends with a summary of results, 
conclusions and proposals for future research. 


2 Model 

We propose a following model of a Call Centre: the analyzed period is finite 
(e.g. one working day) with the system starting empty. The state variable 
X(t) represents total number of service requests (served/waiting calls) in the 
system at time t. The size n{t) of the system, which represents the number of 


an- 
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non empty possible states, is finite, equal to s{t) = number of identical servers 
(agents) plus q{t) = capacity of the queue, with corresponding discrete state 
space ip{t) = {0,n(t)},|(/ 3 (t)| = 1 + s{t) + q{t). Customers arrive according 
to an inhomogenous Poisson process with rate A(t), the service time is i.i.d. 
exponentially distributed with rate /r(t). The load p{t) = X{t)/s{t)p{t) can be 
bigger than 1 . 

Service requests that are not served immediately can leave the system (hang 
up or balk) with probability I- 7 , otherwise, after joining the queue, they aban¬ 
don after reaching their patience time. The patience times are independent and 
identically exponentially distributed with mean l/rj. Queued requests are FCFS 
served. All of this is modeled via the state transition rates of a CTMC which 
is described by infinitesimal generator matrix Q{t) : n(t) -fix n{t) -f 1, Q{i) = 
{qi,j{t)) and the initial state probability vector p( 0 ), where the time dependent 
value qij{i j) is the rate at which the state i changes to the state j and 
qi^i = — Qi,j represents the rate for the event of staying in the same state. 

Because X(t) = k is a birth-and-death process, it can be described by fol¬ 
lowing state dependent birth qk,k+i{t) = Xk{t) and death qk,k-i(t) = fik{t) 
rates: 


Xk{t) 

ti-kit) 


X{t), if 0 < fc < s{t) — 1 
jXft), if s{t) < k < n — 1 

kfi(t), if 1 < fc < s{t) — 1 

s{t)p,{t) + (k — s(t))r], if s{t) < k < n 


( 1 ) 

( 2 ) 


Similar to the model Mi in Deslauriers et al. [ 2007 ] . The transient distri¬ 
bution at time t p(t) for a given time dependent generator matrix Q{t) can be 
calculated using Kolmogorov‘s forward equations: 


pff) = p{f)Q{f) (3) 

where the vector p{t) = [po{t), ...,p„(t)] gives probabilities of the system being 
in any of the states at time t. 

As we do not allow blocking or abandonment due to the overflow of the sys¬ 
tem, the capacity of the queue has to be big enough to be considered practically 
infinite, which is insofar realistic, as the cost of setting practically unlimited 
queue space in the telecommunications equipment is negligible nowadays. The 
system size must, in consequence, ensure that the probability of being in the 
state n (blocking or abandoning service requests) is insignificant compared to 
the required computational precision of the whole model. 


3 Multi-Step Uniformization with Steady-State 
Detection 

The infinitesimal generator matrix Q(t) of an inhomogenous continuous-time 
Markov chain (ICTMC) is time dependent and the process is described by mod¬ 
ified Kolmogorov‘s forward equations (§. 

When the changes in generator matrix Q occur in a discrete way at finite 
points of time and all rates are constant during the intervals between them. 
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we could also replace the analyzed ICTMC with a sequence of homogeneous 
systems computing the state probability vectors for consecutive time periods 
recursively using uniformization as proposed e.g. in Rindos et al. [1995 


Gross and Miller 1984 


or in 


In case of a call center, time dependent changes in Q can occur either dis¬ 
cretely due to the changing number of servers or due to changes in the arrival 
rate. Since the forecast and current traffic data in call center Management 
applications are already aggregated with their average values by an arbitrary 


period (e.g. 5, 15 or SOmin), we will further assume, similarly to Ingolfsson et al. 
2010|, Q(t) being accordingly piecewise constant and refer to such consecutive 


time periods of length A with the coresponding homogenous continuous-time 
Markov chains (HCTMCs) as steps. 

Another approach adopting uniformization for time-inhomogenous CTMCs 
introduced by van Dijk |1992 with subsequent improvements by Van Moorsel 


and Wolter [1998 , Arns et al. [2010 and Andreychenko et al. [2010 could 


be used if continuous arrival rates were available, reducing the error of the 
approximation with the average rates. 

Uniformization or Randomization, known since the publication of Jensen in 
1953 and, therefore, often referenced to as Jensen method, is the method of 
choice for computing transient behavior of CTMCs. Many authors compared 
its performance in different applications with the conclusion that it usually out¬ 
performs known differential equation solvers (e.g. Grassmann [1978 , Reibman 
and Trivedi 1988 , Arns et al. [2010 ). To use uniformization we first define the 
matrix 

^ (4) 


P = I 


a 


which for a > maXi{\qi^i\) is a stochastic matrix, 
uniformization rate. Further, let 


/3(at, k) = e 


-at 


k\ 


The value of a is called 


(5) 


be the probability of a Poisson process with rate a to generate k events in the 
interval [0,t). One now finds for p{t) 


p{t) = p{0) k)(Pf (6) 

fc =0 


The formula (§ can be interpreted as a discrete time Markov process (DTMC) 
embedded in a Poisson process generating events at rate a. 

The implemented uniformization algorithm is based on [Reibman and Trivedi] 


1988 and computes transient state probabilities for a CTMC with the following 
modification of (§ : 

p(t)= Vn(z)e-“‘^ (7) 


2=0 


where a is uniformization rate, as described in Q, and n(i) is the state prob¬ 
ability vector of the underlying DTMC after each step i computed iteratively 
by: 

n(0)=p(0), n(*) = n(*-i)p ( 8 ) 
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To compute p{i), within prespecified error tolerance, in finite time, the compu¬ 
tation stops when the remaining value of cdf of Poisson distribution is less than 
the error bound e: 

i=0 

with k being the right truncation point. As at increases, the corresponding 
probabilities of small number of i Poisson events occurring become less signif¬ 
icant. This allows us to start the summation from the 1‘th. iteration called left 
truncation point with the equation reduced to: 

k 

p(t)=^n(*)e-“‘^ (10) 

i—l 


Reibman and Trivedi |1988| suggests that the values of I and k be derived by: 


E 






i =0 


i\ - 2 


( 11 ) 


The main computational effort of the algorithm lies in consecutive k matrix 
vector multiplications (MVM), necessary for calculation of epochs of DTMC 
in ([^, and is of Olrjk) where rj is the number of nonzero elements of (sparse) 
P. For large at, as the distribution converges to normal, both left and right 


truncation points I and k in (111 will tend to be symmetric to the mean. The 
number 1^ is consequently of 0{at) and the number of additional MVMs 
for the given error tolerance of O'/ai and proportional to inverse cdf for that 
given e. Therefore, although we could solve the p{t) with any accuracy e > 0, 
choosing a higher, acceptable for a respective practical application, value would 
bring some computational advantage. 

The savings due to (tighter) left truncation are, however, rather insignificant, 
unless the computation of the first significant DTMC is performed in a more 
efficient way. 

An example of this, presented first in Muppala and Trivedi [1992 , is based 
on recognizing the steady-state of the underlying DTMC. If convergence of the 
probability vector in (|^ is guaranteed then we can stop the MVM after arriving 
at the steady-state, i.e. let us assume that DTMC has the steady state solution 
n(oo) and that after the S iteration of ([^ ||n( 5 ') — n(oo)||„ < S{S), where ||.||„ 
is an arbitrary vector norm. Then (101 changes to: 


U{S) 

s 


pit) = 


.-at 


n(^)(i-^ 


„-at (o^ty 


i—l 


same as p{t) in (10) 


i=0 


iiS<L 


) if / < S' < A:, (12) 


if S > fc 


with p[t) used instead p(t) denoting transient state probability vector computed 
using approximate steady state DTMC vector n(S). According to Malhotra 
et al. |1994| for a predefined error bound e (as in ([^,([TI|) the following inequality 
holds: 


\\pit)-m\\<l + ms) 


(13) 
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The computing of consecutive epochs of the DTMC is equivalent to the 
power method of finding stationary probability vector of a finite Markov chain. 
According to Stewart |2009 if the stochastic matrix P is aperiodic convergence of 
the power method is guaranteed and the number of iterations k needed to satisfy 
a tolerance criterion ^ may be obtained approximately from the relationship 


P = i-e., k = 


log^ 

logp 


where p is the magnitude of subdominant eigenvalue A 2 of matrix P 
l = ||Ai||>||A2||>||A3||...>|iA^|| 


(14) 


(15) 


reducing, consequently, the computational complexity to 0(ry log^/log\\ 2 \)■ 
Since in most cases the size of the subdominant eigenvalue is not known in 
advance, the usual method of testing for convergence is to examine some norm 
of the difference of successive iterates: 


\\Ili{k) - Ili{k - m)\\ < ^ 


(16) 


Stewart 2009 recommends using the relative convergence test of iterates spaced 


apart by m being function of the rate of convergence: 


maXi 


|ni(fc) - Xli{k - m)\ 

|n.(fc)| 




(17) 


and suggests envisaging a "battery" of different convergence tests in order to 
accept the approximation 11 (S') as being sufficiently accurate. The main risk in 
this approach is that in order to ensure, with the above proposed methods, that 
the n(S) is steady, an additional computational effort for both the convergence 
tests and the required additional number of iterations can easily obliterate the 
potential savings. 

However, in case of our model, we can easily calculate precise stationary dis¬ 
tribution n(oo) in advance, using global balance equations (e.g. as in Stewart 
2009 ) with birth and death rates as in Q and ([^. Therefore, we can conse¬ 


quently, as proposed in Burak |2014| instead of iterating the DTMC vector in 


up to a point S where it would probably satisfy required convergence tests, 
simply use the n((X)) (instead of n(S'), as proposed in the original algorithm by 
Muppala and Trivedi |1992| ) as the f){t -I- A) approximation of pit + A). 

This can be decided after relatively few i iterations due to convergence prop¬ 
erties of the power method as described e.g. in O'Leary et al. [I97^ or in 


standard books on numerical analysis, using numerically estimated convergence 
function of n(i) (as proposed in Burak |2014| ), as it allows for precise calculation 
of the error of such a solution: 


£t-l-A — 


Hn(0-n(oo)||oo 

||n(oo)||oo 


(18) 


in order to decide if it is acceptable (smaller than a predefined steady-detection 
threshold 5t). 

One of the biggest advantages of the uniformization is its strict error bound¬ 
ing for one step independently of its length. It is not difficult to show (e.g. 
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Van Moorsel and Sanders |1997|) that the total error for a number of 


formization steps is the sum of truncation errors (error bounds) for each step. 

Assume for a time period T with a known initial distribution p(0) that for 
any p(r), r = (0, T] the value of each its state has to be computed with an error 
less than Et- Let us further assume Et < Et being the error after computing 
some p{t),t < T. Then: 


£t + cA; < — T — t 


(19) 


As the error bound of steady state approximation is, in case the steady state 
is reached, absolute and independent of the error of the previous steps, we can 
set the convergence threshold dependent rather on the actual total error bound 
than the error for the single step (as proposed e.g. by Malhotra et al. |1994 ). It 
allows, consequently, to trade the error bounds of steps for higher convergence 
thresholds while still within the global error bound for the whole solution. Then, 
assuming the system at time m, 0 < m < T - to satisfy Et < Et for each p{t), 
t = (m, T] we have to: 


< Et — £m — ^ CA 


( 20 ) 


4 Computational Examples 

To test the implementation the following model has been used: a service system 
(call center) working for time T = 24h and starting empty. The arrival rate 
changes sinusoidal with two peaks and is divided into 288 (5min) periods with 
constant averaged rates, same as the first example in [Burak |2014| . The service 
rate and number of servers are constant (/i(t) = p, s{t) = s), the arrival rate 
varies in time - X{t) = s/r(0.85 + 0 . 2 sm( 37 rt/T), 0 < t < T (the load varying 
between 0.65 and 1.05 as shown in Figure]^. The probability 1 —7 of a customer 
immediately leaving when not served immediately is 0.03. The mean value of 
patience time I /77 is equal to 4min. 

The capacity of the queue is constant and chosen so that for all times the 
probability Pn{t) of the system being in the state n is less than 1 x 10“^ for 
all tested system sizes. To evaluate the impact of the proposed steady-state 
detection algorithm, models of 5 different sizes have been at first calculated 
using unmodified uniformization algorithm with an error step e = 1.5 x 10“® 
corresponding to the total error bound et = 2.88 x 10“^. 


Table 1: Computation times, steady-state detection (avx). 


eA=le-7 

S =0(eA=le-5) 

Et ^ 

5e-03 

Et ^ 

1.5e-02 

Et ^ 

3e-02 

Et ^ 

5e-02 

System size 

time 

tjn^ 

time 

tlri^ 

time 


time 

t/n^ 

time 

tjrX' 

54 .(30+24) 

4.25 

1.46 

4.27 

1.46 

3.25 

1.11 

2.13 

0.731 

2.48 

0.850 

150...(100+50) 

15.8 

0.70 

15.3 

0.68 

11.5 

0.51 

4.46 

0.198 

3.73 

0.166 

390... (300+90) 

90.1 

0.59 

86.2 

0.57 

62.0 

0.41 

43.6 

0.286 

15.3 

0.101 

1200(1000+200) 

709 

0.49 

715 

0.50 

534 

0.37 

468 

0.325 

378 

0.262 

3300(3000+300) 

5996 

0.55 

5547 

0.51 

4350 

0.40 

4053 

0.372 

3915 

0.359 


load 0.65 < p < 1.05 


The detailed results of computation times are in Table 
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Figure 1: Computational example - System load 



All experiments were performed on a 1.7GHz PC under 64bit Linux OS with 
a processor supporting vector operations in both: avx with 256bit vectors (4 
double or 8 float operations simultaneously) and the older sse instruction set 
with 128bit vector operations (an Intel i5-3317U with cpu throttling disabled 
via kernel scaling governor), compiled with GNU GCC compiler. 

All measurements use standard Unix time.h/clock() function - returning 
CPU time. All times are in milliseconds. The impact of reduced computa- 


Figure 2: Number of iterations (mvm) per step, system size 1200. 



tional effort due to steady-state detection for some chosen total error bounds 
(between 0 and 5 x 10“^), with corresponding steady-state detection thresholds, 
is illustrated for the system of size 1200 in Figure 

Figure [^shows the expected state of the system, derived from the calculated 
probability vector as: 


ES{t) = '^ip^{t), p{t) = [po-Pn] (21) 

i 

Figure shows its relative error for different steady-state detection thresh¬ 
olds. The reference for the error estimate has been calculated with ca = 
1 X 10-13. 

Figure shows the probability for an incoming service request to be served 
immediately (with no waiting time). 





















































Figure 3: Expected system state, system size 1200. 



Figure 4: Error of the expected system state, system size 1200. 



Figure 5: Probability of a request being served immediately, system size 1200. 



Table shows computation times and maximal values of p„(t),0 < t < T 
- probability of the system being in the state n, for different values of: queue 
length, probability of the customer entering the queue despite not being served 
immediately 7 and mean value of patience time 1 / 77 . The difference compared 
to the corresponding value for 7 = 0.97 and l/? 7 = 4min from the Table (equal 
to 468ms) is not only both due to the bigger size of the system and higher 
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Table 2: Computation time and p(n) depending on queue length, 7 and 1/ry . 


ea = le-7, Et = 3e-02 

7=0.97, 1 / 7 ]= Smin 

7=0.99, 1/?/= Smin 

7=0.997, l/r;=12min 

System size 

time 

max pn{t) 

time 

max pn(t) 

time 

max pn{t) 

1250..(1000+250) 

479 

6.6 X 10-s 

661 

3.8 X 10-’’ 

686 

1.8 X lO-i 

1300.. (1000+300) 

530 

5.8 X 10“i2 

715 

9.3 X 10““ 

760 

7.7 X 10-° 

1400..(1000+400) 

607 

9.2 X 10-2° 

827 

1.1 X 10-1° 

813 

8.7 X 10-1° 


load 0.65 < p < 1.05 


uniformization rate a resulting from higher queue length, but also to some extent 
due to the higher variability of the system state, resulting in fewer steps where 
the steady-state could be detected within the respective threshold. To illustrate 
this effect, we repeated the experiment with the load variability reduced to only 
0.95 < p < 1.05 i.e. with the arrival rate X{t) = s/r(1.0 -|- 0.05sm(37rt/T), 0 < 
t < T. The results corresponding to the cases from the Table are shown in 
the Table m 


Table 3: Computation time and p(n) depending on queue length, 7 and 1/ry . 


ea = le-7, Et = 3e-02 

7=0.97, 1 / 7 ]= 8mm 

7=0.99, 1 / 7 ]= Smin 

7=0.997, l/r;=12min 

System size 

time 

max pn{t) 

time 

max pn(t) 

time 

max p„{t) 

1250..(1000+250) 

34.6 

6.6 X 10-° 

93.6 

3.8 X 10-1 

320 

1.8 X lO-i 

1300..(1000+300) 

34.4 

5.8 X 10-12 

98.8 

9.3 X 10-1° 

349 

7.7 X 10-° 

1400..(1000+400) 

38.0 

9.2 X 10-2° 

no 

1.1 X 10-1° 

380 

8.7 X 10-1° 


load 0.95 < p < 1.05 


5 Conclusion 

In this paper we showed that the uniformization with steady-state detection 
can be used in a very effective way to evaluate transient behavior of multiserver 
queues. Applied to the modeling of the call center schedules, it allows calculation 
of transient system states for systems of any, possible in practical applications, 
size in a very short time, in a numerically stable way, with very high precision, 
using relatively common and inexpensive CPU. It can, therefore, be used for 
schedule planning based on available forecasts, as described in |Ingolfsson et al.| 

2M0] . 

The presented method can be extended in several directions. One could 
be, in regard to call center modeling, to automatically optimize the model size 
(queue length) with significant impact on the computational efficiency. Another 
could be to use known periodicity of traffic forecasts to divide total error bound 
in between known times of the day, bounded by the points of time when the 
system will reach a steady state, than for the whole modeled period. 
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